Предобработка данных

In [1]:
# Подключение к Google drive

from google.colab import drive
drive.mount('/content/drive')
---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
<ipython-input-1-45e8ec640228> in <module>
      1 # Подключение к Google drive
      2 
----> 3 from google.colab import drive
      4 drive.mount('/content/drive')

ModuleNotFoundError: No module named 'google'
In [1]:
import numpy as np
import pandas as pd

Специфические преобразования

In [2]:
# Загрузим набор данных

#df = pd.read_csv('/content/drive/My Drive/Colab Notebooks/freMPL-R.csv', low_memory=False)
df = pd.read_csv('D:\Cloud\Git\geekbrains-ml-business\\freMPL-R.csv', low_memory=False)
df = df.loc[df.Dataset.isin([5, 6, 7, 8, 9])]
df.drop('Dataset', axis=1, inplace=True)
df.dropna(axis=1, how='all', inplace=True)
df.drop_duplicates(inplace=True)
df.reset_index(drop=True, inplace=True)

В предыдущем уроке мы заметили отрицательную величину убытка для некоторых наблюдений. Заметим, что для всех таких полисов переменная "ClaimInd" принимает только значение 0. Поэтому заменим все соответствующие значения "ClaimAmount" нулями.

In [3]:
NegClaimAmount = df.loc[df.ClaimAmount < 0, ['ClaimAmount','ClaimInd']]
print('Unique values of ClaimInd:', NegClaimAmount.ClaimInd.unique())
NegClaimAmount.head()
Unique values of ClaimInd: [0]
Out[3]:
ClaimAmount ClaimInd
82 -74.206042 0
175 -1222.585196 0
177 -316.288822 0
363 -666.758610 0
375 -1201.600604 0
In [4]:
df.loc[df.ClaimAmount < 0, 'ClaimAmount'] = 0

Для моделирования частоты убытков сгенерируем показатель как сумму индикатора того, что убыток произошел ("ClaimInd") и количества заявленных убытков по различным видам ущерба за 4 предшествующих года ("ClaimNbResp", "ClaimNbNonResp", "ClaimNbParking", "ClaimNbFireTheft", "ClaimNbWindscreen").

В случаях, если соответствующая величина убытка равняется нулю, сгенерированную частоту также обнулим.

In [5]:
df['ClaimsCount'] = df.ClaimInd + df.ClaimNbResp + df.ClaimNbNonResp + df.ClaimNbParking + df.ClaimNbFireTheft + df.ClaimNbWindscreen
df.loc[df.ClaimAmount == 0, 'ClaimsCount'] = 0
df.drop(["ClaimNbResp", "ClaimNbNonResp", "ClaimNbParking", "ClaimNbFireTheft", "ClaimNbWindscreen"], axis=1, inplace=True)
In [6]:
pd.DataFrame(df.ClaimsCount.value_counts()).rename({'ClaimsCount': 'Policies'}, axis=1)
Out[6]:
Policies
0.0 104286
2.0 3529
1.0 3339
3.0 2310
4.0 1101
5.0 428
6.0 127
7.0 26
8.0 6
9.0 2
11.0 1
In [7]:
import plotly.express as px
fig = px.scatter(df, x='ClaimsCount', y='ClaimAmount', title='Зависимость между частотой и величиной убытков')
fig.show()

Для моделирования среднего убытка можем рассчитать его как отношение величины убытков к их частоте.

In [8]:
df.loc[df.ClaimsCount > 0, 'AvgClaim'] = df['ClaimAmount']/df['ClaimsCount']

Общие преобразования

Класс для общих случаев

In [9]:
class InsDataFrame:


    ''' Load data method '''

    def load_pd(self, pd_dataframe):
        self._df = pd_dataframe


    ''' Columns match method '''

    def columns_match(self, match_from_to):
        self._df.rename(columns=match_from_to, inplace=True)


    ''' Person data methods '''

    # Gender
    _gender_dict = {'Male':0, 'Female':1}

    def transform_gender(self):
        self._df['Gender'] = self._df['Gender'].map(self._gender_dict)

        

    # Age
    @staticmethod
    def _age(age, age_max):
        if pd.isnull(age):
            age = None
        elif age < 18:
            age = None
        elif age > age_max:
            age = age_max
        return age
      
    def transform_age(self, age_max=70):
        self._df['driver_minage'] = self._df['driver_minage'].apply(self._age, args=(age_max,))

    # Age M/F
    @staticmethod
    def _age_gender(age_gender):
        _age = age_gender[0]
        _gender = age_gender[1]
        if _gender == 0: #Male
            _driver_minage_m = _age
            _driver_minage_f = 18
        elif _gender == 1: #Female
            _driver_minage_m = 18
            _driver_minage_f = _age
        else:
            _driver_minage_m = 18
            _driver_minage_f = 18
        return [_driver_minage_m, _driver_minage_f]
    
    def transform_age_gender(self):
        self._df['driver_minage_m'],self._df['driver_minage_f'] = zip(*self._df[['driver_minage','Gender']].apply(self._age_gender, axis=1).to_frame()[0])

    # Experience
    @staticmethod
    def _exp(exp, exp_max):
        if pd.isnull(exp):
            exp = None
        elif exp < 0:
            exp = None
        elif exp > exp_max:
            exp = exp_max
        return exp

    def transform_exp(self, exp_max=52):
        self._df['driver_minexp'] = self._df['driver_minexp'].apply(self._exp, args=(exp_max,))


    ''' Other data methods '''

    def polynomizer(self, column, n=2):
        if column in list(self._df.columns):
            for i in range(2,n+1):
                self._df[column+'_'+str(i)] = self._df[column]**i

    def get_dummies(self, columns):
        self._df = pd.get_dummies(self._df, columns=columns)


    ''' General methods '''

    def info(self):
        return self._df.info()

    def head(self, columns, n=5):
        return self._df.head(n)

    def len(self):
        return len(self._df)

    def get_pd(self, columns):
        return self._df[columns]

Создаем класс-наследник, в котором переопределяем некоторые методы, специфические для конкретной ситуации, и создаем новые

  • В данных стаж "LicAge" измеряется в неделях.
  • Фактор "SocioCateg" содержит информацию о социальной категории в виде кодов классификации CSP. Агрегируем имеющиеся коды до 1 знака, а затем закодируем их с помощью one-hot encoding.

Wiki

Более подробный классификатор

In [10]:
class InsDataFrame_Fr(InsDataFrame):

    # Experience (weeks to years)
    @staticmethod
    def _exp(exp, exp_max):
        if pd.isnull(exp):
            exp = None
        elif exp < 0:
            exp = None
        else:
            exp * 7 // 365
        if exp > exp_max:
            exp = exp_max
        return exp

    # Marital status
    _MariStat_dict = {'Other':0, 'Alone':1}

    def transform_MariStat(self):
        self._df['MariStat'] = self._df['MariStat'].map(self._MariStat_dict)
    
    # Social category
    def transform_SocioCateg(self):
        self._df['SocioCateg'] = self._df['SocioCateg'].str.slice(0,4)

Подгружаем данные

In [11]:
data = InsDataFrame_Fr()
In [12]:
data.load_pd(df)
In [13]:
data.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 115155 entries, 0 to 115154
Data columns (total 17 columns):
Exposure       115155 non-null float64
LicAge         115155 non-null int64
RecordBeg      115155 non-null object
RecordEnd      59455 non-null object
Gender         115155 non-null object
MariStat       115155 non-null object
SocioCateg     115155 non-null object
VehUsage       115155 non-null object
DrivAge        115155 non-null int64
HasKmLimit     115155 non-null int64
BonusMalus     115155 non-null int64
ClaimAmount    115155 non-null float64
ClaimInd       115155 non-null int64
OutUseNb       115155 non-null float64
RiskArea       115155 non-null float64
ClaimsCount    115155 non-null float64
AvgClaim       10869 non-null float64
dtypes: float64(6), int64(5), object(6)
memory usage: 14.9+ MB

Преобразовываем параметры

In [14]:
# Переименовываем
data.columns_match({'DrivAge':'driver_minage','LicAge':'driver_minexp'})
In [15]:
# Преобразовываем
data.transform_age()
data.transform_exp()
data.transform_gender()
data.transform_MariStat()
data.transform_SocioCateg()
In [16]:
# Пересечение пола и возраста, их квадраты
data.transform_age_gender()
data.polynomizer('driver_minage_m')
data.polynomizer('driver_minage_f')

Для переменных, содержащих более 2 значений:

  • либо упорядочиваем значения,
  • либо используем фиктивные переменные (one-hot encoding).

NB: В H2O не рекомендуется использовать one-hot encoding. Тем не менее, используем здесь фиктивные переменные, чтобы в дальнейшем сохранить возможность сравнения результатов построенных моделей.

In [17]:
# Onehot encoding
data.get_dummies(['VehUsage','SocioCateg'])
In [18]:
data.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 115155 entries, 0 to 115154
Data columns (total 30 columns):
Exposure                           115155 non-null float64
driver_minexp                      115155 non-null int64
RecordBeg                          115155 non-null object
RecordEnd                          59455 non-null object
Gender                             115155 non-null int64
MariStat                           115155 non-null int64
driver_minage                      115155 non-null int64
HasKmLimit                         115155 non-null int64
BonusMalus                         115155 non-null int64
ClaimAmount                        115155 non-null float64
ClaimInd                           115155 non-null int64
OutUseNb                           115155 non-null float64
RiskArea                           115155 non-null float64
ClaimsCount                        115155 non-null float64
AvgClaim                           10869 non-null float64
driver_minage_m                    115155 non-null int64
driver_minage_f                    115155 non-null int64
driver_minage_m_2                  115155 non-null int64
driver_minage_f_2                  115155 non-null int64
VehUsage_Private                   115155 non-null uint8
VehUsage_Private+trip to office    115155 non-null uint8
VehUsage_Professional              115155 non-null uint8
VehUsage_Professional run          115155 non-null uint8
SocioCateg_CSP1                    115155 non-null uint8
SocioCateg_CSP2                    115155 non-null uint8
SocioCateg_CSP3                    115155 non-null uint8
SocioCateg_CSP4                    115155 non-null uint8
SocioCateg_CSP5                    115155 non-null uint8
SocioCateg_CSP6                    115155 non-null uint8
SocioCateg_CSP7                    115155 non-null uint8
dtypes: float64(6), int64(11), object(2), uint8(11)
memory usage: 17.9+ MB
In [19]:
col_features = [
                'driver_minexp',
                'Gender',
                'MariStat',
                'HasKmLimit',
                'BonusMalus',
                'OutUseNb',
                'RiskArea',
                'driver_minage_m',
                'driver_minage_f',
                'driver_minage_m_2',
                'driver_minage_f_2',
                'VehUsage_Private',
                'VehUsage_Private+trip to office',
                'VehUsage_Professional',
                'VehUsage_Professional run',
                'SocioCateg_CSP1',
                'SocioCateg_CSP2',
                'SocioCateg_CSP3',
                'SocioCateg_CSP4',
                'SocioCateg_CSP5',
                'SocioCateg_CSP6',
                'SocioCateg_CSP7'  
]
In [20]:
col_target = ['ClaimAmount', 'ClaimsCount', 'AvgClaim']
In [21]:
df_freq = data.get_pd(col_features+col_target)
In [22]:
df_ac = df_freq[df_freq['ClaimsCount'] > 0].reset_index().copy()

Разделение набора данных на обучающую, валидационную и тестовую выборки

In [23]:
from sklearn.model_selection import train_test_split
In [24]:
# Разбиение датасета для частоты на train/val/test

x_train_c, x_test_c, y_train_c, y_test_c = train_test_split(df_freq[col_features], df_freq.ClaimsCount, test_size=0.3, random_state=1)
x_valid_c, x_test_c, y_valid_c, y_test_c = train_test_split(x_test_c, y_test_c, test_size=0.5, random_state=1)
In [25]:
# Разбиение датасета для среднего убытка на train/val/test 

x_train_ac, x_test_ac, y_train_ac, y_test_ac = train_test_split(df_ac[col_features], df_ac.AvgClaim, test_size=0.3, random_state=1)
x_valid_ac, x_test_ac, y_valid_ac, y_test_ac = train_test_split(x_test_ac, y_test_ac, test_size=0.5, random_state=1)

Обобщенные линейные модели (Generalized Linear Models, GLM)

Теория

Пусть $y$ – целевая переменная, $X$ – матрица объясняющих переменных, $\beta$ – вектор параметров модели.

Матрица $X$ составлена из всех векторов наблюдений $x_i$, каждый из которых представляет собой объясняющую переменную.

Основные компоненты обобщенной линейной модели:

  • Систематическая компонента $\eta$:
    • $\eta = X\beta \hspace{10pt}(=\beta_0+\beta_1x_1+\beta_2x_2+\dots+\beta_nx_n)$.
  • Случайная компонента $y$:
    • Элементы вектора $y$ – независимые одинаково распределенные случайные величины, имеющие функцию плотности распределения $f(y;\theta,\phi)$ из экспоненциального семейства.
    • Распределения из экспоненциального семейства имеют параметры $\theta$ (характеристика среднего) и $\phi$ (характеристика дисперсии). В общем виде данные распределения могут быть определены: $$f_i(y_i;\theta_i,\phi)=\exp\left\lbrace \frac{y_i\theta_i-b(\theta_i)}{a_i(\phi)} + c(y_i, \phi) \right\rbrace,$$ где $a_i(\phi)$, $b(\theta_i)$ и $c(y_i, \phi)$ некоторые функции.
    • Для распределений из данного семейства дисперсия является функцией от среднего.
    • Экспоненциальное семейство включает распределения нормальное, экспоненциальное, Пуассона, гамма, хи-квадрат, бета и другие.
  • Функция связи $g$:
    • $\mathbb{E}\left[y\right]=\mu=g^{-1}\left(\eta\right)$, $\mu$ – математическое ожидание $y$;
    • $g$ – монотонная дифференцируемая функция.

GLM с распределением Пуассона

Регрессия Пуассона обычно используется в случаях, когда зависимая переменная представляет собой счетные значения и ошибки предполагаются распределенными в соответствии с распределением Пуассона. Зависимая переменная должна быть неотрицательной.

Функции связи между таргетом и объясняющими переменными предполагается логарифмической: $$g(\eta) = \ln(\eta) \Rightarrow \hat {y} = e^{x{^T}\beta + {\beta_{0}}}.$$

Модель оценивается методом максимального правдоподобия, функция логарифма правдоподобия с учетом штрафа регуляризации эластичной сети имеет вид: $$\max_{\beta,\beta_0} \frac{1}{N} \sum_{i=1}^{N} \Big( y_i(x_{i}^{T}\beta + \beta_0) - e^{x{^T_i}\beta + {\beta_0}} \Big)- \lambda \Big( \alpha||\beta||_1 + \dfrac {1} {2}(1 - \alpha)||\beta||^2_2 \Big),$$

где

  • $\lambda$ – параметр, отвечающий за силу регуляризации. $\lambda\in\mathbb{R}^{+}$;
  • $\alpha$ – параметр, отвечающий за распределение штрафов регуляризации между нормой 1 ($\ell_1$) и нормой 2 ($\ell_2$). $\alpha\in[0,1]$;
  • $||\beta||{_1}$ – штраф регуляризации $\ell_1$ (LASSO). $||\beta||{_1} = \sum{^p_{k=1}} |\beta{_k}|$;
  • $||\beta||{_2}$ – штраф регуляризации $\ell_2$ (Ridge). $||\beta||{_2} = \sum{^p_{k=1}} \beta{^2_k}$.

Тогда соответствующая метрика Deviance имеет вид: $$D = -2 \sum_{i=1}^{N} \big( y_i \text{ln}(y_i / \hat {y}_i) - (y_i - \hat {y}_i) \big).$$

Вывод функции правдоподобия для GLM с распределением Пуассона (без регуляризации)

NB: Ниже $\lambda$ не имеет отношения к вышеупомянутому одноименному параметру регуляризации.

Напомним, что функция вероятности для распределения Пуассона имеет вид: $$p(k;\lambda) = \frac{\lambda^k e^{-\lambda}}{k!}, \hspace{10pt} \lambda\in\mathbb{R}^{+}.$$ Также, для распределения Пуассона справедливо, что: $$\mathbb{E}\left[k\right] = Var(k) = \lambda.$$ Тогда, для оценивания коэффициентов нашей модели необходимо максимизировать правдоподобие (совместную условную вероятность при имеющихся данных), что данные имеют распределение Пуассона: $$p(y_1,\dots,y_n|x_1,\dots,x_n;\beta_0,\beta) = \prod_{i=1}^{N}\frac{e^{y_i(x_i{^T}\beta + {\beta_{0}})} e^{-e^{x_i{^T}\beta + {\beta_{0}}}}}{y_i!} = L(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n).$$ Для упрощения задачи оптимизации перейдем к логарифму правдоподобия: $$\ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n) = \sum_{i=1}^{N}\left(y_i(x_i{^T}\beta + {\beta_{0}}) -e^{x_i{^T}\beta + {\beta_{0}}}-\ln(y_i!)\right).$$ Поскольку величина $\ln(y_i!)$ не зависит от выбора параметров, можно упростить задачу: $$\ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n) = \sum_{i=1}^{N}\left(y_i(x_i{^T}\beta + {\beta_{0}}) -e^{x_i{^T}\beta + {\beta_{0}}}\right).$$ Далее решается задача оптимизации для определения параметров модели (пакеты реализуют оптимизацию численными методами): $$\frac{\partial \ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n)}{\partial \beta_0} = 0,\\\frac{\partial \ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n)}{\partial \beta} = 0.$$ Обычно минимизируется отрицательное правдоподобие, которое является выпуклой функцией.

Вывод метрики Deviance для GLM с распределением Пуассона

Метрика Deviance представляет собой отношение правдоподобия между двумя моделями: рассматриваемой моделью и "идеальной" моделью, в которая бы идеально предсказывала бы зависимую переменную. $$Deviance = 2(\ell_{ideal} - \ell_{model})$$

В качестве такой "идеальной модели" может использоваться сама зависимая переменная. Тогда, логарифм правдоподобия "идеальной модели" для GLM с распределением Пуассона имеет вид: $$\ell_{ideal} = \sum_{i=1}^{N}\left(y_i \ln(y_i) -y_i-\ln(y_i!)\right).$$

Из приведенного выше вывода правдоподобия для рассматриваемой модели, мы можем записать, обозначив $\hat{y}_i = e^{x{^T}\beta + {\beta_{0}}}$: $$\ell_{model} = \sum_{i=1}^{N}\left(y_i \ln(\hat{y}_i) -\hat{y}_i-\ln(y_i!)\right).$$

Тогда получаем, $$Deviance = 2\sum_{i=1}^{N}\left(y_i \ln(y_i) -y_i - y_i \ln(\hat{y}_i) +\hat{y}_i\right) = -2\sum_{i=1}^{N}\left(y_i \ln(y_i/\hat{y}_i) - (y_i -\hat{y}_i)\right).$$

GLM с гамма-распределением

GLM с гамма-распределением используется для моделирования положительной непрерывной зависимой переменной, когда ее условная дисперсия увеличивается вместе со средним значением, но коэффициент вариации зависимой переменной предполагается постоянным.

Обычно GLM с гамма-распределением используются с логарифмической или обратной функциями связи: $$g(\eta) = \ln(\eta);\hspace{20pt}g(\eta) = \frac{1}{\eta}.$$

Модель оценивается методом максимального правдоподобия, функция логарифма правдоподобия (для обратной функции связи) с учетом штрафа регуляризации эластичной сети имеет вид: $$\max_{\beta,\beta_0} - \frac{1}{N} \sum_{i=1}^{N} \frac{y_i}{x{^T_i}\beta + \beta_0} + \text{ln} \big( x{^T_i}\beta + \beta_0 \big ) - \lambda \Big( \alpha||\beta||_1 + \dfrac {1} {2}(1 - \alpha)||\beta||^2_2 \Big),$$

где

  • $\lambda$ – параметр, отвечающий за силу регуляризации. $\lambda\in\mathbb{R}^{+}$;
  • $\alpha$ – параметр, отвечающий за распределение штрафов регуляризации между нормой 1 ($\ell_1$) и нормой 2 ($\ell_2$). $\alpha\in[0,1]$;
  • $||\beta||{_1}$ – штраф регуляризации $\ell_1$ (LASSO). $||\beta||{_1} = \sum{^p_{k=1}} |\beta{_k}|$;
  • $||\beta||{_2}$ – штраф регуляризации $\ell_2$ (Ridge). $||\beta||{_2} = \sum{^p_{k=1}} \beta{^2_k}$.

Соответствующая метрика Deviance имеет вид: $$D = 2 \sum_{i=1}^{N} - \text{ln} \bigg (\dfrac {y_i} {\hat {y}_i} \bigg) + \dfrac {(y_i - \hat{y}_i)} {\hat {y}_i}.$$

Вывод функции правдоподобия для GLM с Гамма-распределением и логарифмической функцией связи (без регуляризации)

NB: Ниже $\alpha$ и $\lambda$ не имеют отношения к вышеупомянутым одноименным параметрам регуляризации.

Напомним, что функция плотности вероятности для Гамма-распределения имеет вид: $$f(x;\alpha,\lambda) = \begin{cases} \frac{\alpha^\lambda}{\Gamma(\lambda)}x^{\lambda-1}e^{-\alpha x},&x \ge 0\\ 0,& x <0 \end{cases},\hspace{20pt} \Gamma(x) = \int_0^{\infty}x^{\lambda-1}e^{-x}dx. $$ Также, для Гамма-распределения справедливо, что: $$ \mathbb{E}[x] = \frac{\lambda}{\alpha},\hspace{15pt} Var(x) = \frac{\lambda}{\alpha^2}.$$

Для оценивания GLM удобно параметризовать данное распределение иначе: $$ \mathbb{E}[x] = \frac{\lambda}{\alpha} = \mu,\hspace{15pt} Var(x) = \frac{\lambda}{\alpha^2} = \frac{\mu^2}{\lambda}.$$ Тогда, $$f(x;\mu,\lambda) = \frac{1}{\Gamma(\lambda)}\alpha^\lambda x^{\lambda-1}e^{-\alpha x} = \frac{1}{\Gamma(\lambda)}\left(\frac{\lambda}{\mu}\right)^\lambda x^{\lambda-1}e^{-\left(\frac{\lambda}{\mu}\right)x} = \frac{1}{x\cdot\Gamma(\lambda)}\left(\frac{\lambda x}{\mu}\right)^\lambda e^{-\frac{\lambda x}{\mu}},\hspace{15pt} x\ge 0.$$

Также напомним, что для GLM с функцией связи $\ln(x)$ предсказание (оценка математического ожидания) имеет вид $\hat {y} = e^{x{^T}\beta + {\beta_{0}}}$.

Тогда, для оценивания коэффициентов нашей модели необходимо максимизировать правдоподобие (совместную условную вероятность при имеющихся данных), что данные имеют Гамма-распределение, в предположении известного параметра $\lambda$: $$p(y_1,\dots,y_n|x_1,\dots,x_n;\beta_0,\beta) = \prod_{i=1}^{N}\frac{1}{y_i\cdot\Gamma(\lambda)}\left(\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right)^\lambda e^{-\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}} = L(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n).$$ Для упрощения задачи оптимизации перейдем к логарифму правдоподобия: \begin{align} \ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n) &= \sum_{i=1}^{N}\left(-\ln(\Gamma(\lambda))-\ln(y_i)+ \lambda\ln\left(\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right) -\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right)\\ &= \sum_{i=1}^{N}\left(-\ln(\Gamma(\lambda))-\ln(y_i) + \lambda\ln\left(\lambda\right)+\lambda\ln\left(y_i\right)-\lambda\left({x_i{^T}\beta + {\beta_{0}}}\right) -\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right)\\ &= \sum_{i=1}^{N}\left((\lambda-1)\ln\left(y_i\right)-\ln(\Gamma(\lambda)) + \lambda\ln\left(\lambda\right)-\lambda\left({x_i{^T}\beta + {\beta_{0}}}\right) -\frac{\lambda y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right). \end{align} Тогда обозначим слагаемые, не зависящие от параметров модели ($\beta_0$, $\beta$), как некоторую константу $C$, получаем: $$\ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n) = \sum_{i=1}^{N}\left(-\lambda\left(\left({x_i{^T}\beta + {\beta_{0}}}\right) +\frac{y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right) + C(y_i, \lambda)\right),$$ Таким образом, требуется максимизировать $$\ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n) = -\sum_{i=1}^{N}\left(\left({x_i{^T}\beta + {\beta_{0}}}\right) + \frac{y_i}{e^{x_i{^T}\beta + {\beta_{0}}}}\right),$$ Далее решается задача оптимизации для определения параметров модели (пакеты реализуют оптимизацию численными методами): $$\frac{\partial \ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n)}{\partial \beta_0} = 0,\\\frac{\partial \ell(\beta_0,\beta|y_1,\dots,y_n,x_1,\dots,x_n)}{\partial \beta} = 0.$$ Обычно минимизируется отрицательное правдоподобие, которое является выпуклой функцией.

Вывод метрики Deviance для GLM с Гамма-распределением

Метрика Deviance представляет собой отношение правдоподобия между двумя моделями: рассматриваемой моделью и "идеальной" моделью, в которая бы идеально предсказывала бы зависимую переменную. $$Deviance = 2(\ell_{ideal} - \ell_{model})$$

В качестве такой "идеальной модели" может использоваться сама зависимая переменная. Тогда, логарифм правдоподобия "идеальной модели" для GLM с распределением Пуассона имеет вид: $$\ell_{ideal} = -\sum_{i=1}^{N}\left(1 + \ln(y_i)\right).$$

Из приведенного выше вывода правдоподобия для рассматриваемой модели, мы можем записать, обозначив $\hat{y}_i = e^{x{^T}\beta + {\beta_{0}}}$: $$\ell_{model} = -\sum_{i=1}^{N}\left(\frac{y_i}{\hat{y}_i} + \ln(\hat{y}_i)\right).$$

Тогда получаем, $$Deviance = 2\sum_{i=1}^{N}\left(- 1 - \ln(y_i) + \frac{y_i}{\hat{y}_i} + \ln(\hat{y}_i)\right) = 2\sum_{i=1}^{N}\left(-\ln\left(\frac{y_i}{\hat{y}_i}\right)+\frac{y_i-\hat{y}_i}{\hat{y}_i} \right).$$

Дополнительная литература по GLM

Установка H2O на Google Colaboratory и инициализация

In [26]:
!apt-get install default-jre
"apt-get" ­Ґ пў«пҐвбп ў­гв७­Ґ© Ё«Ё ў­Ґи­Ґ©
Є®¬ ­¤®©, ЁбЇ®«­пҐ¬®© Їа®Ја ¬¬®© Ё«Ё Ї ЄҐв­л¬ д ©«®¬.
In [27]:
!java -version
"java" ­Ґ пў«пҐвбп ў­гв७­Ґ© Ё«Ё ў­Ґи­Ґ©
Є®¬ ­¤®©, ЁбЇ®«­пҐ¬®© Їа®Ја ¬¬®© Ё«Ё Ї ЄҐв­л¬ д ©«®¬.
In [28]:
!pip install h2o
^C
In [26]:
import h2o
from h2o.estimators.glm import H2OGeneralizedLinearEstimator
h2o.init()
Checking whether there is an H2O instance running at http://localhost:54321 ..... not found.
Attempting to start a local H2O server...
; Java HotSpot(TM) 64-Bit Server VM (build 25.202-b08, mixed mode)
  Starting server from C:\Users\mmingalov\AppData\Local\Continuum\anaconda3\lib\site-packages\h2o\backend\bin\h2o.jar
  Ice root: C:\Users\MMINGA~1\AppData\Local\Temp\tmpurc959hk
  JVM stdout: C:\Users\MMINGA~1\AppData\Local\Temp\tmpurc959hk\h2o_mmingalov_started_from_python.out
  JVM stderr: C:\Users\MMINGA~1\AppData\Local\Temp\tmpurc959hk\h2o_mmingalov_started_from_python.err
  Server is running at http://127.0.0.1:54321
Connecting to H2O server at http://127.0.0.1:54321 ... successful.
H2O_cluster_uptime: 03 secs
H2O_cluster_timezone: Asia/Vladivostok
H2O_data_parsing_timezone: UTC
H2O_cluster_version: 3.30.1.1
H2O_cluster_version_age: 6 days
H2O_cluster_name: H2O_from_python_mmingalov_sm5m9l
H2O_cluster_total_nodes: 1
H2O_cluster_free_memory: 3.526 Gb
H2O_cluster_total_cores: 8
H2O_cluster_allowed_cores: 8
H2O_cluster_status: accepting new members, healthy
H2O_connection_url: http://127.0.0.1:54321
H2O_connection_proxy: {"http": null, "https": null}
H2O_internal_security: False
H2O_API_Extensions: Amazon S3, Algos, AutoML, Core V3, TargetEncoder, Core V4
Python_version: 3.7.4 final

Построение GLM для частоты страховых случаев

In [27]:
# Преобразование в H2O-Frame

h2o_train_c = h2o.H2OFrame(pd.concat([x_train_c, y_train_c], axis=1))
h2o_valid_c = h2o.H2OFrame(pd.concat([x_valid_c, y_valid_c], axis=1))
h2o_test_c = h2o.H2OFrame(pd.concat([x_test_c, y_test_c], axis=1))
Parse progress: |█████████████████████████████████████████████████████████| 100%
Parse progress: |█████████████████████████████████████████████████████████| 100%
Parse progress: |█████████████████████████████████████████████████████████| 100%
In [28]:
# Инициализируем и обучим GLM модель c кросс-валидацией

glm_poisson = H2OGeneralizedLinearEstimator(family = "poisson", link = "Log", nfolds=5)
glm_poisson.train(y="ClaimsCount", x = h2o_train_c.names[1:-1], training_frame = h2o_train_c, validation_frame = h2o_valid_c)
glm Model Build progress: |███████████████████████████████████████████████| 100%
In [29]:
# Параметры модели: распределение, функция связи, гиперпараметры регуляризации, количество использованных объясняющих переменных

glm_poisson.summary()
GLM Model: summary
family link regularization number_of_predictors_total number_of_active_predictors number_of_iterations training_frame
0 poisson log Elastic Net (alpha = 0.5, lambda = 9.667E-5 ) 21 19 3 Key_Frame__upload_9e5dedfd95f071da184ad208e5594569.hex
Out[29]:

In [30]:
# Метрики качества модели - по всем данным и на кросс-валидации

glm_poisson.cross_validation_metrics_summary().as_data_frame()
Out[30]:
mean sd cv_1_valid cv_2_valid cv_3_valid cv_4_valid cv_5_valid
0 mae 0.39096567 0.0014827132 0.39165637 0.39039108 0.39292508 0.3909263 0.38892958
1 mean_residual_deviance 1.0610672 0.011107087 1.0652161 1.0592242 1.0725753 1.0652282 1.043092
2 mse 0.5885074 0.007385119 0.5875963 0.5841486 0.6003421 0.58951 0.58094
3 null_deviance 17457.518 158.76501 17517.955 17314.377 17617.734 17573.697 17263.826
4 r2 0.008276904 0.0020501784 0.008408169 0.008374793 0.008159138 0.005325864 0.011116553
5 residual_deviance 17106.238 218.06065 17149.979 16963.477 17277.043 17331.264 16809.428
6 rmse 0.7671303 0.004802894 0.76654834 0.7642961 0.77481747 0.76779556 0.7621942
7 rmsle 0.35898024 0.0018329391 0.35989913 0.3587962 0.36080462 0.35942647 0.35597476
In [31]:
# Таблица коэффициентов модели (в зависимости от модели могут выводиться также стандартная ошибка, z-score и p-value)

glm_poisson._model_json['output']['coefficients_table'].as_data_frame()
Out[31]:
names coefficients standardized_coefficients
0 Intercept -2.565659 -1.575867
1 Gender -0.163708 -0.079349
2 MariStat -0.114494 -0.041245
3 HasKmLimit -0.439625 -0.137290
4 BonusMalus 0.012823 0.196838
5 OutUseNb 0.081984 0.057037
6 RiskArea 0.010119 0.022421
7 driver_minage_m 0.000209 0.003928
8 driver_minage_f 0.012369 0.198390
9 driver_minage_m_2 0.000000 0.000000
10 driver_minage_f_2 -0.000091 -0.113901
11 VehUsage_Private -0.190581 -0.090123
12 VehUsage_Private+trip to office 0.000000 0.000000
13 VehUsage_Professional 0.320065 0.105780
14 VehUsage_Professional run 0.491840 0.065987
15 SocioCateg_CSP1 -0.079219 -0.012003
16 SocioCateg_CSP2 -0.205329 -0.034207
17 SocioCateg_CSP3 0.242698 0.025177
18 SocioCateg_CSP4 0.050770 0.012696
19 SocioCateg_CSP5 -0.003721 -0.001769
20 SocioCateg_CSP6 0.060435 0.024822
21 SocioCateg_CSP7 -0.495631 -0.005237
In [32]:
# Таблица нормированных коэффициентов по всем данным и на кросс-валидации

pmodels = {}
pmodels['overall'] = glm_poisson.coef_norm()
for x in range(len(glm_poisson.cross_validation_models())):
    pmodels[x] = glm_poisson.cross_validation_models()[x].coef_norm()
coef = pd.DataFrame.from_dict(pmodels).round(5)
coef['overall'] = abs(coef['overall'])
coef.sort_values('overall',ascending=False)
Out[32]:
overall 0 1 2 3 4
Intercept 1.57587 -1.57828 -1.57496 -1.58069 -1.57913 -1.56726
driver_minage_f 0.19839 0.15420 0.23662 0.12252 0.08085 0.08563
BonusMalus 0.19684 0.19391 0.19892 0.19069 0.20767 0.18914
HasKmLimit 0.13729 -0.13715 -0.14122 -0.14004 -0.14488 -0.12421
driver_minage_f_2 0.11390 -0.06694 -0.14951 -0.05228 -0.01813 -0.02039
VehUsage_Professional 0.10578 0.10345 0.10427 0.11369 0.09605 0.09829
VehUsage_Private 0.09012 -0.08717 -0.09123 -0.08745 -0.11763 -0.08690
Gender 0.07935 -0.07461 -0.08685 -0.05681 -0.04077 -0.06523
VehUsage_Professional run 0.06599 0.06905 0.06056 0.06666 0.06598 0.06189
OutUseNb 0.05704 0.06167 0.05393 0.06353 0.05544 0.05082
MariStat 0.04125 -0.03616 -0.02935 -0.03896 -0.04968 -0.05435
SocioCateg_CSP2 0.03421 -0.03344 -0.04498 -0.02751 -0.03846 -0.02947
SocioCateg_CSP3 0.02518 0.02561 0.01731 0.02435 0.03259 0.02597
SocioCateg_CSP6 0.02482 0.01823 0.02953 0.02185 0.03669 0.01081
RiskArea 0.02242 0.01754 0.03416 0.01859 0.02606 0.01562
SocioCateg_CSP4 0.01270 0.02055 0.00000 0.01493 0.00949 0.01765
SocioCateg_CSP1 0.01200 -0.01601 -0.00566 -0.01451 -0.01367 -0.00979
SocioCateg_CSP7 0.00524 -0.00455 -0.00405 -0.02964 -0.00262 -0.00310
driver_minage_m 0.00393 0.01478 0.02019 0.00261 0.01214 0.00000
SocioCateg_CSP5 0.00177 -0.00080 -0.00638 0.00000 -0.00068 0.00000
driver_minage_m_2 0.00000 0.00000 -0.01602 0.00000 0.00237 -0.00862
VehUsage_Private+trip to office 0.00000 0.00000 0.00000 0.00000 -0.02034 0.00000
In [33]:
# Построение прогнозных значений для обучающей, валидационной и тестовой выборок

c_train_pred = glm_poisson.predict(h2o_train_c).as_data_frame()
c_valid_pred = glm_poisson.predict(h2o_valid_c).as_data_frame()
c_test_pred = glm_poisson.predict(h2o_test_c).as_data_frame()
glm prediction progress: |████████████████████████████████████████████████| 100%
glm prediction progress: |████████████████████████████████████████████████| 100%
glm prediction progress: |████████████████████████████████████████████████| 100%
In [34]:
# Сохранение обученной модели

model_glm_poisson = h2o.save_model(model=glm_poisson, path="/content/drive/My Drive/Colab Notebooks/", force=True)
In [35]:
model_glm_poisson
Out[35]:
'C:\\content\\drive\\My Drive\\Colab Notebooks\\GLM_model_python_1597634047802_1'

Построение GLM для среднего убытка

In [36]:
# Преобразование в H2O-Frame

h2o_train_ac = h2o.H2OFrame(pd.concat([x_train_ac, y_train_ac], axis=1))
h2o_valid_ac = h2o.H2OFrame(pd.concat([x_valid_ac, y_valid_ac], axis=1))
h2o_test_ac = h2o.H2OFrame(pd.concat([x_test_ac, y_test_ac], axis=1))
Parse progress: |█████████████████████████████████████████████████████████| 100%
Parse progress: |█████████████████████████████████████████████████████████| 100%
Parse progress: |█████████████████████████████████████████████████████████| 100%
In [37]:
# Инициализируем и обучим GLM модель c кросс-валидацией

glm_gamma = H2OGeneralizedLinearEstimator(family = "gamma", link = "Log", nfolds=5)
glm_gamma.train(y="AvgClaim", x = h2o_train_ac.names[1:-1], training_frame = h2o_train_ac, validation_frame = h2o_valid_ac)
glm Model Build progress: |███████████████████████████████████████████████| 100%
In [48]:
# Параметры модели: распределение, функция связи, гиперпараметры регуляризации, количество использованных объясняющих переменных

glm_gamma.summary()
GLM Model: summary
family link regularization number_of_predictors_total number_of_active_predictors number_of_iterations training_frame
0 gamma log Elastic Net (alpha = 0.5, lambda = 2.602E-4 ) 21 19 3 Key_Frame__upload_9fd96e4015fbd6e7da40e3ce77624407.hex
Out[48]:

In [49]:
# Метрики качества модели - по всем данным и на кросс-валидации

glm_gamma.cross_validation_metrics_summary().as_data_frame()
Out[49]:
mean sd cv_1_valid cv_2_valid cv_3_valid cv_4_valid cv_5_valid
0 mae 1190.5558 51.509735 1162.4398 1176.5247 1149.2976 1185.1139 1279.4027
1 mean_residual_deviance 2.097358 0.12192993 1.9642259 2.0791116 2.1277795 2.029054 2.2866187
2 mse 1.552649E7 2.036847E7 9533856.0 6720484.5 4624020.0 4957996.5 5.1796092E7
3 null_deviance 3178.4329 211.89409 3044.8904 3238.793 3074.522 3010.4248 3523.5347
4 r2 -0.0033354152 0.006333948 0.0053234217 -0.0062541133 -0.011550041 -0.003892568 -3.0377542E-4
5 residual_deviance 3191.71 214.81517 3003.3013 3280.8381 3119.3247 3031.4065 3523.6794
6 rmse 3450.8088 2126.7373 3087.6943 2592.3896 2150.3535 2226.656 7196.95
7 rmsle 1.7673236 0.050552726 1.731823 1.753846 1.8565875 1.7465256 1.7478359
In [50]:
# Таблица коэффициентов модели (в зависимости от модели могут выводиться также стандартная ошибка, z-score и p-value)

glm_gamma._model_json['output']['coefficients_table'].as_data_frame()
Out[50]:
names coefficients standardized_coefficients
0 Intercept 7.310145 7.017288
1 Gender -0.314623 -0.152506
2 MariStat 0.299338 0.110533
3 HasKmLimit -0.066859 -0.016827
4 BonusMalus -0.002345 -0.039579
5 OutUseNb -0.021085 -0.015749
6 RiskArea 0.030396 0.067861
7 driver_minage_m -0.021733 -0.397791
8 driver_minage_f 0.012161 0.187589
9 driver_minage_m_2 0.000179 0.268589
10 driver_minage_f_2 -0.000154 -0.183893
11 VehUsage_Private 0.000000 0.000000
12 VehUsage_Private+trip to office 0.024089 0.012019
13 VehUsage_Professional -0.021973 -0.008023
14 VehUsage_Professional run -0.290266 -0.044476
15 SocioCateg_CSP1 -0.129253 -0.020325
16 SocioCateg_CSP2 -0.218260 -0.037997
17 SocioCateg_CSP3 -0.069526 -0.006869
18 SocioCateg_CSP4 -0.019099 -0.005290
19 SocioCateg_CSP5 0.000000 0.000000
20 SocioCateg_CSP6 0.092551 0.035955
21 SocioCateg_CSP7 -1.358740 -0.015578
In [51]:
# Таблица нормированных коэффициентов по всем данным и на кросс-валидации

pmodels = {}
pmodels['overall'] = glm_gamma.coef_norm()
for x in range(len(glm_gamma.cross_validation_models())):
    pmodels[x] = glm_gamma.cross_validation_models()[x].coef_norm()
coef_ac = pd.DataFrame.from_dict(pmodels).round(5)
coef_ac['overall'] = abs(coef['overall'])
coef_ac.sort_values('overall',ascending=False)
Out[51]:
overall 0 1 2 3 4
Intercept 1.57587 7.02409 7.01198 7.03483 7.00690 6.99880
driver_minage_f 0.19839 0.33528 0.00412 0.22712 0.40513 0.00000
BonusMalus 0.19684 -0.01186 -0.06135 -0.03188 -0.06618 -0.02150
HasKmLimit 0.13729 -0.01963 -0.02375 -0.01289 -0.02306 -0.00253
driver_minage_f_2 0.11390 -0.30378 0.00000 -0.19363 -0.39472 -0.05700
VehUsage_Professional 0.10578 -0.00069 -0.03386 -0.00457 -0.02320 0.00704
VehUsage_Private 0.09012 0.01080 0.00000 0.00000 0.00146 -0.02314
Gender 0.07935 -0.13295 -0.16333 -0.18032 -0.12736 -0.16659
VehUsage_Professional run 0.06599 -0.03698 -0.04013 -0.08061 -0.04393 -0.03116
OutUseNb 0.05704 -0.00505 -0.01090 -0.02326 -0.01646 -0.01679
MariStat 0.04125 0.11567 0.13478 0.12144 0.13968 0.03105
SocioCateg_CSP2 0.03421 -0.02897 -0.02917 -0.05004 -0.04984 -0.03261
SocioCateg_CSP3 0.02518 -0.01939 0.00000 0.00000 -0.00965 -0.00629
SocioCateg_CSP6 0.02482 0.07343 0.00366 0.01611 0.05233 0.02250
RiskArea 0.02242 0.05435 0.07500 0.06091 0.05421 0.09358
SocioCateg_CSP4 0.01270 -0.00126 -0.02119 -0.01324 0.00256 0.00258
SocioCateg_CSP1 0.01200 -0.02192 -0.01640 -0.01849 -0.01259 -0.03095
SocioCateg_CSP7 0.00524 -0.01805 -0.01839 -0.01406 -0.01753 0.00000
driver_minage_m 0.00393 -0.14504 -0.47552 -0.42862 -0.09956 -0.87043
SocioCateg_CSP5 0.00177 0.00000 0.00000 0.00000 0.00000 -0.00936
driver_minage_m_2 0.00000 0.04616 0.34372 0.31696 0.00000 0.66458
VehUsage_Private+trip to office 0.00000 0.00000 0.00478 0.03661 0.00000 0.00000
In [52]:
# Построение прогнозных значений для обучающей, валидационной и тестовой выборок

ac_train_pred = glm_gamma.predict(h2o_train_ac).as_data_frame()
ac_valid_pred = glm_gamma.predict(h2o_valid_ac).as_data_frame()
ac_test_pred = glm_gamma.predict(h2o_test_ac).as_data_frame()
glm prediction progress: |████████████████████████████████████████████████| 100%
glm prediction progress: |████████████████████████████████████████████████| 100%
glm prediction progress: |████████████████████████████████████████████████| 100%
In [53]:
# Сохранение обученной модели

model_glm_gamma = h2o.save_model(model=glm_gamma, path="/content/drive/My Drive/Colab Notebooks/", force=True)
In [54]:
model_glm_gamma
Out[54]:
'C:\\content\\drive\\My Drive\\Colab Notebooks\\GLM_model_python_1597634047802_2'

Использование GLM моделей

In [55]:
h2o_df = h2o.H2OFrame(df_freq[col_features])
Parse progress: |█████████████████████████████████████████████████████████| 100%
In [56]:
df['CountPredicted'] = glm_poisson.predict(h2o_df).as_data_frame()
df['AvgClaimPredicted'] = glm_gamma.predict(h2o_df).as_data_frame()
glm prediction progress: |████████████████████████████████████████████████| 100%
glm prediction progress: |████████████████████████████████████████████████| 100%
In [57]:
df['BurningCost'] = df.CountPredicted * df.AvgClaimPredicted
df[['CountPredicted', 'AvgClaimPredicted', 'BurningCost']].head()
Out[57]:
CountPredicted AvgClaimPredicted BurningCost
0 0.268689 1087.936414 292.316704
1 0.268689 1087.936414 292.316704
2 0.240122 1146.063987 275.195142
3 0.200097 1127.364931 225.582287
4 0.168869 1106.883916 186.918109

* Домашнее задание: GLM для прогнозирования наступления страхового случая

In [64]:
# Разбиение датасета на train/val/test

df_ind = data.get_pd(col_features+['ClaimInd'])

x_train_ind, x_test_ind, y_train_ind, y_test_ind = train_test_split(df_ind[col_features], df_ind.ClaimInd, test_size=0.3, random_state=1)
x_valid_ind, x_test_ind, y_valid_ind, y_test_ind = train_test_split(x_test_ind, y_test_ind, test_size=0.5, random_state=1)
In [65]:
# Преобразование в H2O-Frame
h2o_train_ind = h2o.H2OFrame(pd.concat([x_train_ind, y_train_ind], axis=1))
h2o_valid_ind = h2o.H2OFrame(pd.concat([x_valid_ind, y_valid_ind], axis=1))
h2o_test_ind = h2o.H2OFrame(pd.concat([x_test_ind, y_test_ind], axis=1))
h2o_train_ind
Parse progress: |█████████████████████████████████████████████████████████| 100%
Parse progress: |█████████████████████████████████████████████████████████| 100%
Parse progress: |█████████████████████████████████████████████████████████| 100%
driver_minexp Gender MariStat HasKmLimit BonusMalus OutUseNb RiskArea driver_minage_m driver_minage_f driver_minage_m_2 driver_minage_f_2 VehUsage_Private VehUsage_Private+trip to office VehUsage_Professional VehUsage_Professional run SocioCateg_CSP1 SocioCateg_CSP2 SocioCateg_CSP3 SocioCateg_CSP4 SocioCateg_CSP5 SocioCateg_CSP6 SocioCateg_CSP7 ClaimInd
52 1 0 0 50 0 10 18 55 324 3025 0 1 0 0 0 0 0 0 1 0 0 1
52 0 0 0 50 0 2 70 18 4900 324 1 0 0 0 0 0 0 0 0 1 0 1
52 1 1 0 90 1 10 18 61 324 3721 0 1 0 0 0 0 0 0 1 0 0 0
52 1 0 0 80 0 10 18 26 324 676 0 1 0 0 0 0 0 0 1 0 0 0
52 0 0 0 50 1 10 70 18 4900 324 1 0 0 0 0 0 0 0 0 1 0 0
52 0 0 0 72 0 4 30 18 900 324 0 1 0 0 0 0 0 0 1 0 0 0
52 0 0 0 54 0 9 42 18 1764 324 0 1 0 0 0 0 0 0 1 0 0 0
52 0 0 0 50 0 4 70 18 4900 324 1 0 0 0 0 0 0 0 0 1 0 0
52 0 0 0 50 0 6 34 18 1156 324 0 1 0 0 0 0 0 1 0 0 0 0
52 0 0 0 72 1 11 48 18 2304 324 0 1 0 0 0 0 0 1 0 0 0 0
Out[65]:

In [66]:
# Преобразуем целевую переменную ClaimInd в категориальную при помощи метода asfactor во всех наборах данных

h2o_train_ind['ClaimInd'] = h2o_train_ind['ClaimInd'].asfactor()
h2o_valid_ind['ClaimInd'] = h2o_valid_ind['ClaimInd'].asfactor()
h2o_test_ind['ClaimInd'] = h2o_test_ind['ClaimInd'].asfactor()
h2o_train_ind
driver_minexp Gender MariStat HasKmLimit BonusMalus OutUseNb RiskArea driver_minage_m driver_minage_f driver_minage_m_2 driver_minage_f_2 VehUsage_Private VehUsage_Private+trip to office VehUsage_Professional VehUsage_Professional run SocioCateg_CSP1 SocioCateg_CSP2 SocioCateg_CSP3 SocioCateg_CSP4 SocioCateg_CSP5 SocioCateg_CSP6 SocioCateg_CSP7 ClaimInd
52 1 0 0 50 0 10 18 55 324 3025 0 1 0 0 0 0 0 0 1 0 0 1
52 0 0 0 50 0 2 70 18 4900 324 1 0 0 0 0 0 0 0 0 1 0 1
52 1 1 0 90 1 10 18 61 324 3721 0 1 0 0 0 0 0 0 1 0 0 0
52 1 0 0 80 0 10 18 26 324 676 0 1 0 0 0 0 0 0 1 0 0 0
52 0 0 0 50 1 10 70 18 4900 324 1 0 0 0 0 0 0 0 0 1 0 0
52 0 0 0 72 0 4 30 18 900 324 0 1 0 0 0 0 0 0 1 0 0 0
52 0 0 0 54 0 9 42 18 1764 324 0 1 0 0 0 0 0 0 1 0 0 0
52 0 0 0 50 0 4 70 18 4900 324 1 0 0 0 0 0 0 0 0 1 0 0
52 0 0 0 50 0 6 34 18 1156 324 0 1 0 0 0 0 0 1 0 0 0 0
52 0 0 0 72 1 11 48 18 2304 324 0 1 0 0 0 0 0 1 0 0 0 0
Out[66]:

In [67]:
# Инициализируем и обучим GLM модель c кросс-валидацией
glm_binomial = H2OGeneralizedLinearEstimator(family = "binomial", link = "logit", nfolds=5)
glm_binomial.train(y="ClaimInd", x = h2o_train_ind.names[1:-1], training_frame = h2o_train_ind, validation_frame = h2o_valid_ind)
glm Model Build progress: |███████████████████████████████████████████████| 100%
In [68]:
# Параметры модели: распределение, функция связи, гиперпараметры регуляризации, количество использованных объясняющих переменных
glm_binomial.summary()
GLM Model: summary
family link regularization number_of_predictors_total number_of_active_predictors number_of_iterations training_frame
0 binomial logit Elastic Net (alpha = 0.5, lambda = 2.368E-5 ) 21 21 3 py_4_sid_8387
Out[68]:

In [69]:
# Метрики качества модели - по всем данным и на кросс-валидации
glm_binomial.cross_validation_metrics_summary().as_data_frame()
Out[69]:
mean sd cv_1_valid cv_2_valid cv_3_valid cv_4_valid cv_5_valid
0 accuracy 0.48834112 0.123260565 0.5701125 0.3071154 0.5769897 0.41211557 0.57537234
1 auc 0.5632411 0.006997295 0.5660236 0.5548154 0.5567837 0.56824255 0.5703405
2 aucpr 0.11385259 0.003487019 0.11366126 0.114658564 0.1085848 0.11401406 0.118344285
3 err 0.5116589 0.123260565 0.42988747 0.69288456 0.42301026 0.5878844 0.42462766
4 err_count 8248.6 1987.1151 6953.0 11179.0 6835.0 9462.0 6814.0
5 f0point5 0.13276389 0.0038562226 0.13376741 0.12888747 0.13326208 0.12945792 0.1384446
6 f1 0.18737845 0.0025702408 0.18573603 0.18845735 0.18446486 0.1871134 0.19112061
7 f2 0.3199156 0.02258622 0.30373833 0.35041305 0.29956597 0.3373606 0.30850005
8 lift_top_group 1.4511671 0.29972917 1.6584636 1.2513475 1.0589594 1.799441 1.487624
9 logloss 0.31126416 0.0051384293 0.30746928 0.31919405 0.30856302 0.3073845 0.31370994
10 max_per_class_error 0.56132007 0.12132138 0.4730897 0.74869436 0.4870604 0.6201439 0.47761193
11 mcc 0.058239557 0.0054304483 0.059452318 0.049687397 0.056770682 0.063738786 0.061548598
12 mean_per_class_accuracy 0.5479004 0.006927682 0.5507276 0.535893 0.5482588 0.55292803 0.5516945
13 mean_per_class_error 0.4520996 0.006927682 0.44927236 0.46410698 0.45174125 0.44707194 0.4483055
14 mse 0.08540736 0.0018239913 0.084031336 0.08817648 0.08434692 0.08410823 0.08637381
15 null_deviance 10102.831 155.4547 10013.807 10354.365 10019.579 9975.458 10150.945
16 pr_auc 0.11385259 0.003487019 0.11366126 0.114658564 0.1085848 0.11401406 0.118344285
17 precision 0.11119933 0.0043023136 0.11273813 0.10645452 0.11245272 0.107396446 0.11695482
18 r2 0.0039267405 0.0010766101 0.0042756135 0.0029711605 0.0026125074 0.0047661806 0.0050082407
19 recall 0.6217437 0.142363 0.5269103 0.8204804 0.51293963 0.726 0.52238804
20 residual_deviance 10036.041 160.32864 9946.017 10299.753 9971.522 9894.707 10068.207
21 rmse 0.29223213 0.003109846 0.2898816 0.29694524 0.29042542 0.2900142 0.29389423
22 specificity 0.4740571 0.15167534 0.57454497 0.25130567 0.58357793 0.3798561 0.581001
In [70]:
# Таблица коэффициентов модели (в зависимости от модели могут выводиться также стандартная ошибка, z-score и p-value)
glm_binomial._model_json['output']['coefficients_table'].as_data_frame()
Out[70]:
names coefficients standardized_coefficients
0 Intercept -2.475348e+00 -2.279640
1 Gender -6.416985e-02 -0.031103
2 MariStat -6.469029e-02 -0.023304
3 HasKmLimit -3.698560e-01 -0.115502
4 BonusMalus 6.792484e-03 0.104267
5 OutUseNb 6.145252e-02 0.042753
6 RiskArea 9.240660e-03 0.020474
7 driver_minage_m -3.727295e-03 -0.070187
8 driver_minage_f -1.168549e-03 -0.018743
9 driver_minage_m_2 -6.145260e-06 -0.009595
10 driver_minage_f_2 -9.652992e-08 -0.000121
11 VehUsage_Private -1.437545e-01 -0.067979
12 VehUsage_Private+trip to office -1.656392e-02 -0.008276
13 VehUsage_Professional 2.260091e-01 0.074695
14 VehUsage_Professional run 2.965853e-01 0.039791
15 SocioCateg_CSP1 -1.073627e-01 -0.016267
16 SocioCateg_CSP2 -1.462515e-01 -0.024365
17 SocioCateg_CSP3 1.043504e-01 0.010825
18 SocioCateg_CSP4 2.983183e-02 0.007460
19 SocioCateg_CSP5 -2.803307e-02 -0.013326
20 SocioCateg_CSP6 4.853385e-02 0.019934
21 SocioCateg_CSP7 -2.029539e-01 -0.002144
In [71]:
# Таблица нормированных коэффициентов по всем данным и на кросс-валидации
pmodels = {}
pmodels['overall'] = glm_binomial.coef_norm()
for x in range(len(glm_binomial.cross_validation_models())):
    pmodels[x] = glm_binomial.cross_validation_models()[x].coef_norm()
coef_ac = pd.DataFrame.from_dict(pmodels).round(5)
coef_ac['overall'] = abs(coef['overall'])
coef_ac.sort_values('overall',ascending=False)
Out[71]:
overall 0 1 2 3 4
Intercept 1.57587 -2.27489 -2.29090 -2.27687 -2.27425 -2.28277
driver_minage_f 0.19839 -0.00833 -0.02781 -0.01141 -0.06330 -0.01514
BonusMalus 0.19684 0.10182 0.09818 0.11236 0.09340 0.10610
HasKmLimit 0.13729 -0.12061 -0.11459 -0.11988 -0.10672 -0.11551
driver_minage_f_2 0.11390 0.00000 0.00748 -0.01283 0.02720 0.00475
VehUsage_Professional 0.10578 0.08324 0.09200 0.06905 0.08850 0.06898
VehUsage_Private 0.09012 -0.04630 -0.06198 -0.06887 -0.05738 -0.07284
Gender 0.07935 -0.06001 -0.06970 -0.03197 -0.05169 -0.02606
VehUsage_Professional run 0.06599 0.04148 0.05271 0.02995 0.05047 0.03569
OutUseNb 0.05704 0.04016 0.04794 0.04378 0.03897 0.04348
MariStat 0.04125 -0.01169 -0.01866 -0.04210 -0.02287 -0.02621
SocioCateg_CSP2 0.03421 -0.02508 -0.01591 -0.01870 -0.03159 -0.01138
SocioCateg_CSP3 0.02518 0.00823 0.01375 0.02643 0.01357 0.00805
SocioCateg_CSP6 0.02482 0.00000 0.02984 0.05772 0.02819 0.02213
RiskArea 0.02242 0.02594 0.02616 0.03645 0.00603 0.00873
SocioCateg_CSP4 0.01270 0.01249 0.00000 0.03382 0.01554 0.00773
SocioCateg_CSP1 0.01200 -0.01363 -0.00854 -0.01187 -0.02037 -0.00829
SocioCateg_CSP7 0.00524 -0.02876 -0.00129 0.00075 0.00024 -0.00082
driver_minage_m 0.00393 -0.23007 -0.24732 -0.04554 -0.26134 -0.03051
SocioCateg_CSP5 0.00177 -0.01415 0.00000 0.02060 0.00000 -0.01165
driver_minage_m_2 0.00000 0.14078 0.13259 -0.03846 0.14751 -0.03201
VehUsage_Private+trip to office 0.00000 0.00000 0.00000 0.00000 0.00000 -0.00439
In [72]:
# Построение прогнозных значений для обучающей, валидационной и тестовой выборок
ind_train_pred = glm_binomial.predict(h2o_train_ind).as_data_frame()
ind_valid_pred = glm_binomial.predict(h2o_valid_ind).as_data_frame()
ind_test_pred = glm_binomial.predict(h2o_test_ind).as_data_frame()
glm prediction progress: |████████████████████████████████████████████████| 100%
glm prediction progress: |████████████████████████████████████████████████| 100%
glm prediction progress: |████████████████████████████████████████████████| 100%
In [75]:
from sklearn.metrics import accuracy_score, f1_score, confusion_matrix
In [76]:
# Выведем импортированные выше метрики классификации для обучающей, валидационной и тестовой выборок
ind_train_pred.head()
Out[76]:
predict p0 p1
0 0 0.907688 0.092312
1 0 0.925845 0.074155
2 1 0.883330 0.116670
3 1 0.885729 0.114271
4 0 0.915998 0.084002
In [ ]:
 
In [77]:
# Выведем импортированные выше метрики классификации для обучающей, валидационной и тестовой выборок
y_train_pred = ind_train_pred['predict']
y_valid_pred = ind_valid_pred['predict']
y_test_pred = ind_test_pred['predict']

print(f'train_accuracy_score: {accuracy_score(y_train_ind, y_train_pred)}')
print(f'valid_accuracy_score: {accuracy_score(y_valid_ind, y_valid_pred)}')
print(f'test_accuracy_score: {accuracy_score(y_test_ind, y_test_pred)}')
train_accuracy_score: 0.579446208813021
valid_accuracy_score: 0.5827013257685405
test_accuracy_score: 0.5784415885145305
In [78]:
print(f'train_f1_score: {f1_score(y_train_ind, y_train_pred)}')
print(f'valid_f1_score: {f1_score(y_valid_ind, y_valid_pred)}')
print(f'test_f1_score: {f1_score(y_test_ind, y_test_pred)}')
train_f1_score: 0.18681635002878524
valid_f1_score: 0.18627229622939714
test_f1_score: 0.18454647256438972
In [79]:
print(f'train_confusion_matrix: {confusion_matrix(y_train_ind, y_train_pred)}')
print(f'valid_confusion_matrix: {confusion_matrix(y_valid_ind, y_valid_pred)}')
print(f'test_confusion_matrix: {confusion_matrix(y_test_ind, y_test_pred)}')
train_confusion_matrix: [[42814 30159]
 [ 3741  3894]]
valid_confusion_matrix: [[9240 6424]
 [ 784  825]]
test_confusion_matrix: [[9168 6481]
 [ 801  824]]

Какие проблемы вы здесь видите? Как можно улучшить данный результат?

Исходя из метрик, модель получилась со слабой предсказательной способностью (и низким значением f-метрики). Возможные причины:

  • мы не провели анализ, отбор и построение признаков по отношению к новой целевой переменной.

Можно улучшить результат с помощью проведения анализа влияния признаков на переменную, а также построения новых признаков с последующим отбором лучших из них.

In [ ]: